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FOREWORD 

This  is  the  Final  Scientific  Report  prepared  by  the  Pratt  & Whitney  Aircraft  Group 
of  United  Technologies  Corporation.  The  effort  was  sponsored  by  the  Air  Force 
Office  of  Scientific  Research,  (AFSC),  United  States  Air  Force,  under  Contract 
F44620-74-C-0060,  project  number  2307-Bl,  with  Mr.  W.  J.  Walker  (NA)  as  Pro- 
ject Engineer.  The  period  covered  is  April  1,  1974  through  September  30,  1977. 

Appreciation  is  expressed  for  the  contributions  made  by  Dr.  L.  H.  Seitelman  and 
Mr.  D.  W.  Snow  of  Pratt  & Whitney  Aircraft  during  the  final  phase  of  the  contract. 
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SECTION  1 
INTRODUCTION 

1.1  SUMMARY 

The  overall  objective  of  this  program  has  been  the  development  of  an  analytical  crack  model- 
ing and  stress  analysis  capability,  based  on  the  boundary-integral  equation  method,  of  suffi- 
cient accuracy  and  efficiency  to  favorably  influence  the  low  cycle  fatigue  life  prediction  pro-  | 

cess  for  gas  turbine  engine  structures.  This  report  discusses  progress  toward  this  objective  un-  \ 

der  the  present  contract,  with  particular  attention  to  the  period  since  the  publication  of  the  ; 

last  Interim  Scientific  Report  in  May  1976.  ; 

The  boundary-integral  equation  (BIE)  method  is  known  to  be  particularly  well  suited  for  the 
calculation  of  the  rapidly  varying  stress  fields  associated  with  both  fracture  mechanics  and 
stress  concentration  problems  (Reference  1).  During  this  program,  the  three-dimensional 
BIE  method  has  been  applied  to  a variety  of  fracture  mechanics  problems.  The  results  of  this 
work  have  been  incorporated  in  the  low  cycle  fatigue  life  prediction  process.  A new  boun- 
dary-integral equation  was  derived  for  use  in  three-dimensional  crack  problems,  which  should 
significantly  improve  the  ability  to  evaluate  such  finite  geometry  effects  as  the  proximity  of 
free  surfaces  or  a condition  of  high  local  surface  curvature.  Most  recently,  modified  BIE  crack 
tip  elements  have  been  developed  and  applied  to  fracture  mechanics  problems. 

The  capabilities  of  the  BIE  method  for  stress  analysis  of  uncracked  gas  turbine  engine  struc- 
tures have  been  expanded  by  preliminary  development  of  a technique  for  merging  BIE  and 
finite  element  analyses,  by  extension  of  the  BIE  method  to  deal  with  anisotropic  and  inho- 
mogeneous materials,  and  by  incorporation  of  thermal  stresses  in  the  BIE  formulation. 

Section  1.2  of  this  report  is  a brief  overview  of  the  BIE  method;  Section  1.3  discusses  the  ap- 
plication of  the  BIE  method  in  the  gas  turbine  engine  environment.  Section  2 is  addressed 
specifically  to  the  question  of  elastic  fracture  mechanics  modeling,  and  Section  3 discusses 
application  of  the  BIE  method  to  inhomogeneous  materials. 

In  addition  to  its  impact  on  Pratt  & Whitney  Aircraft  structural  analysis  efforts,  the  research 
carried  out  under  this  contract  has  been  made  widely  available  to  the  technical  community 
through  contract  reports  (References  2,  3,  and  4),  presentations  at  technical  meetings,  and 
publication  in  the  open  literature.  References  5 through  9 discuss  work  fully  or  partially  sup- 
ported by  this  contract.  In  addition,  the  Appendix  to  this  report  has  been  accepted  for  pub- 
lication in  the  International  Journal  for  Numerical  Methods  in  Engineering.  | 

1.2  OVERVIEW  OF  THE  BOUNDARY-INTEGRAL  EQUATION  METHOD 

The  BIE  method  for  the  solution  of  elastic  stress  analysis  problems  is  based  on  classical  re-  i 

suits  in  mathematical  analysis  and  continuum  mechanics.  It  has  become  practical  as  a solu-  ; 

tion  technique  for  problems  with  general  geometry  and  loading  because  of  the  speed  and 
storage  capacity  of  present  computers. 
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The  analytical  basis  of  the  BIE  method  is  the  replacement,  employing  Betti’s  reciprocal  work 
theorem,  of  the  governing  partial  differential  equations  by  an  integral  identity  for  the  elastic 
displacements. 


uj(p) 


(p,Q)  Uj  (Q)  dS 


+ 


/Uij  (P.O) 
S 


tj  (Q)  dS 


(1) 


In  Eq.  (1),  uj(p)  is  the  displacement  vector  at  an  interior  point  p(x);  tj(Q)  and  Uj(0)  are  the 
boundary  values  of  traction  and  displacement.  The  kernel  functions  Tjj(p,Q)  and  Uj|(p,Q) 
are  the  tractions  and  displacements  in  the  Xj  directions  at  Q(x)  due  to  orthogonal  unit  loads 
in  the  xj  directions  at  p(x). 

Equation  (1)  contains  the  totality  of  boundary  data;  in  general,  a well-posed  elasticity  prob- 
lem consists  of  specifying  values  of  t;  on  part  of  the  boundary,  Sj  and  Uj  on  the  remainder  of 
the  boundary  S^.  Allowing  p(x)  ->■  P(x),  a boundary  point,  Eq.  ( 1 ) becomes  a set  of  integral 
constraint  equations  relating  boundary  displacements  to  boundary  tractions;  following  Refer- 
ences 1 0 and  1 1 , the  BIE  is  obtained. 


Ui(P)/2  + 


(P,Q)  Uj(Q)  dS 


Ujj(P,Q)  tj(Q)  dS 


(2) 


After  the  boundary  solution  uj,  tj  has  been  obtained,  interior  displacements  can  be  calculated 
directly  from  Eq.  (1 ).  Interior  stress  can  be  evaluated  by  differentiation  of  Eq.  (1 ) with 
respect  to  p(x)  and  application  of  Hooke’s  Law.  A detailed  presentation  of  the  analytical 
basis  of  the  BIE  method,  including  an  extensive  bibliography,  was  recently  prepared  (Refer- 
ence 2)  under  support  of  this  contract. 


Solution  of  the  analytical  BIE,  Eq.  (2)  is  possible  for  only  a few  simple  geometries  and  load- 
ings. The  practical  application  of  the  method  is  based  on  two  approximations: 

1.  The  boundary  S of  the  region  to  be  analyzed  is  represented  by  a finite  number  of 
surface  patches,  Sj..  Examples  are  triangles  in  three  dimensions  and  straight  line 
segments  in  two  dimensions. 


2.  The  boundary  data  (uj,  tj)  are  taken  to  have  a known  form  of  variation  on  each 
Sj^;  for  example,  linear  variation  over  each  individual  S|^. 

The  first  assumption  allows  the  BIE  to  be  rewritten  as 


N 

^ifP)  \ 

2 ^ ^ 
k = 1 


TjjfP.Q)  Uj(Q)  dS 


E 

' = 1 


UjjfP.Q)  tj(Q)  dS 


(3) 
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Depending  on  the  particular  variation  chosen  for  geometry  variation  and  boundary  data,  the 
integrals  over  the  surface  patches  are  carried  out  either  numerically  (Reference  12)  or  in 
closed  form  (Reference  13)  leading  to  a linear  algebraic  form  of  the  BIE 

MW  • [T]  |u 

Since  the  BIE  in  both  its  analytical  and  algebraic  forms  involves  only  boundary  data,  no  in- 
terior idealization  is  required.  The  coefficient  matrices  in  Eq.  (4)  are  thus  of  much  smaller 
order  than  those  for  a finite  element  analysis  of  the  same  problem,  although  they  lack  the 
banded  structure  of  finite  element  coefficient  matrices. 

1.3  APPLICATION  OF  BIE  STRESS  ANALYSIS  IN  THE  GAS  TURBINE  ENGINE 
ENVIRONMENT 

BIE  applications  in  the  gas  turbine  engine  environment  fall  into  two  major  classes,  fracture 
mechanics  analysis  (discussed  in  detail  in  Section  2)  and  calculation  of  elastic  stress  fields 
near  notches  and  other  structural  details.  This  section  is  primarily  addressed  to  the  latter 
class,  although  many  of  the  comments  made  also  apply  to  the  fracture  mechanics  problem. 

1.3.1  Geometrical  Complexity 

The  outstanding  feature  of  problems  encountered  in  gas  turbine  engine  applications  is  geome- 
trical complexity,  leading  to  very  large  problem  size.  For  example,  the  modeling  of  turbine 
disk  rims  involves  such  features  as  the  load  bearing  teeth  on  the  disk  lug  and  the  doubly 
curved  surface  at  the  rim  slot/cooling  hole  intersection  (Figure  la). 

The  ability  to  model  parts  using  the  BIE  method  is  strongly  influenced  by  the  way  in  which 
geometry  and  boundary  data  are  allowed  to  vary  over  each  surface  patch.  The  BINTEQ  pro- 
gram, used  in  much  of  the  work  carried  out  under  this  contract,  uses  plane  triangular  surface 
patches  with  linear  boundary  data  variation  (Reference  13).  Modeling  of  a geometry  such  as 
a disk  rim  slot  with  BINTEQ  requires  use  of  a very  large  number  of  elements  to  obtain  ade- 
quate definition  (Figure  lb),  leads  to  long  computer  run  times,  and  makes  the  analysis  un- 
suited to  routine  use. 

An  alternative  (References  1 2 and  14)  is  the  use  of  isoparametric  shape  functions  for  the  re- 
presentation of  both  geometry  and  boundary  data.  The  basis  of  the  isoparametric  method  is 
the  mapping  of  a planar  curve  or  a surface  patch  to  a standard  interval  or  square  by  means  of 
a fixed  set  of  shape  functions.  A particular  set  of  quadratic  shape  functions  is  shown  in  Fig- 
ure 2.  In  this  case  a surface  patch  defined  by  eight  nodes  (x“  , a = 1 , ...,  8)  is  mapped  to  the 
square  (- 1 < { | < 1,-1  < ^2  ^ ^ using  the  relation, 

Xj  (5)  = M“  (t)  Xj"  (5) 

The  variation  of  Xjtf)  is  quadratic  along  the  edges  of  the  .square,  although  cubic  terms  can 
occur  in  the  interior.  Similar  shape  functions  exist  (Reference  14)  for  mapping  of  the  tri- 
angle-like surface  patches  required  for  efficient  mesh  size  transition.  The  improved  modeling 
efficiency  of  the  quadratic  shape  function  approach  (demonstrated  for  the  rim  slot/cooling 
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hole  intersection  in  Figure  lb)  allows  the  same  level  of  solution  accuracy  to  be  achieved 
with  a much  less  refined  surface  map.  Other  familities  of  shape  functions  (linear  and  cubic) 
are  discussed  in  Reference  12. 

Boundary  data  can  be  modeled  using  the  same  representation  as  used  for  the  geometry,  the 
approach  of  Reference  14,  or  different  representations  can  be  allowed  (Reference  12).  The 
con  bination  of  linear  boundary  data  modeling  with  quadratic  geometry  representation  is 
particularly  attractive  since  it  permits  a low-cost,  preliminary  analysis  of  a complex  part 
without  sacrificing  geometric  definition.  If  the  results  of  the  preliminary  analysis  are  satis- 
factory, then  a more  refined  analysis  can  be  carried  out  by  improving  the  boundary  data  re- 
presentation without  changing  the  geometric  representation.  A different  approach  to  the 
geometrical  modeling  problem  is  found  in  substructuring.  Substructuring  refers  to  the  BIE 
modeling  of  a part  in  two  or  more  separate  subregions  which  are  then  joined  by  enforcing 
appropriate  continuity  conditions  on  the  common  elements  which  form  the  interfaces  be- 
tween subregions  (Reference  12). 

Pratt  & Whitney  Aircraft  has  acquired,  for  in-house  use,  a computer  program  (Reference  12), 
using  the  shape  function  approach  and  incorporating  substructuring  capability,  referred  to  as 
BASQUE  (Boundary  Solution  using  Quadratic  Elements).*  This  program  has  been  used  for 
much  of  the  work  carried  out  in  the  current  part  of  the  contract  effort.  Experience  with  the 
BASQUE  code  has  shown  that  the  increase  in  modeling  efficiency  due  to  the  isoparametric 
approach  and  the  substructuring  capability  permits  practical  three-dimensional  BIE  analysis 
of  turbine  disk  rim  and  turbine  blade  attachment  structures. 

Finally,  work  has  been  carried  out  under  this  contract  to  explore  the  technique  of  hybridi- 
zation, that  is,  the  merging  of  BIE  and  finite  element  analysis.  Some  rather  promising  nume- 
rical results  were  achieved  in  the  context  of  two-dimensional  analysis.  These  results  tegether 
with  a discussion  of  more  basic  issues  involved  in  coupling  the  BIE  and  finite  element  meth- 
ods were  reported  in  Reference  3 and  5. 

i.3.2  Loading  Complexity 

Gas  turbine  engine  structures  experience  loads  other  than  the  surface  mechanical  loads  ac- 
counted for  in  Eq.  (1).  In  particular,  the  structures  are  subjected  to  thermal,  and  often  cen- 
trifugal loads.  In  order  to  account  for  these  effects,  the  BIE  must  be  modified  to  include 
the  volume  integrals  of  the  body  force  terms  for  thermal  and  centrifugal  loading.  Since  one 
of  the  advantages  of  the  BIE  method  is  the  absence  of  an  interior  idealization,  the  existence 
of  these  volume  integrals  poses  a difficult  question. 

During  an  earlier  phase  of  this  contract,  a technique  was  developed  for  reducing  the  centrifu- 
gal body  force  volume  integral  to  a surface  integral.  An  exact  reduction  of  the  thermoelastic 
case  to  a surface  integral  was  also  obtained  under  the  assumption  of  a steady  state  tempera- 
ture field.  Both  analyses  were  incorporated  in  BINTEQ  and  verified  numerically  (Reference 
4). 


♦This  proprietary  code  was  developed  at  the  Centre  Technique  des  Industries  Mbchaniques 
(CETIM)  under  the  direction  of  Dr.  J.  C.  Lachat.  Further  information  about  the  code  is 
available  from  M.  Lange,  Department  Calcul  des  Structures,  CETIM,  BP  67,  60404  Senlis, 
France. 
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The  surface  integral  reduction  of  the  thermal  loading  was  also  studied  as  an  approximate  tech- 
nique in  the  case  of  time  dependent  temperature  fields  but  was  found  to  be  unsuitable  for  even 
small  departures  from  steady  state  conditions.  This  is  of  considerable  importance  since,  typi- 
cally, engine  structures  must  be  analyzed  at  several  points  during  an  engine  operating  cycle, 
and  the  temperature  fields  at  most  or  all  of  these  points  result  from  nonsteady  state  condi- 
tions. An  alternative  approach  for  BIE  analysis  of  parts  subjected  to  nonsteady  state  thermal 
loads  is  discussed  in  Section  3. 

1.3.3  Material  Complexity 

Two  types  c!' complexity  arise  in  the  specification  of  elastic  material  properties  for  gas  tur- 
bine engine  structures:  anisotropy  and  inhomogeneity. 

Anisotropy  results  from  the  search  for  new  materials  which  will  have  longer  service  lives 
under  the  increasingly  demanding  operating  conditions  found  in  advanced  engines.  It  has 
been  found  that  the  use  of  anisotropic  materials  ^particularly  directionally  solidified  and 
single  crystal  materials  in  turbine  airfoils)  allows  optimization  of  preferred  material  direc- 
tions with  respect  to  part  loadings. 

Anisotropy  poses  no  problem  in  the  formulation  of  Eq.  (2),  and  integral  representations  for 
the  required  three-dimensional  point  load  functions  were  established  in  References  15  and 
16;  however,  efficient  numerical  evaluation  of  these  point  load  functions  posed  a consider- 
able problem.  During  the  present  contract,  a technique  was  developed  for  minimizing  this 
numerical  problem.  The  technique  was  incorporated  in  the  BASQUE  code  for  verification 
and  study  of  computational  efficiency.  A detailed  description  of  this  work  is  contained  in 
the  Appendix  of  this  report. 

Inhomogeneity  of  elastic  material  properties  in  gas  turbine  engine  parts  nonnally  arises  from 
the  combination  of  temperature-dependent  material  properties  and  the  fact  that  parts  are 
subjected  to  nonuniform  temperature  fields.  This  problem  is  closely  tied  to  the  nonsteady 
state  thermal  body  force  problem  mentioned  above.  One  aspect  of  the  present  contract  ac- 
tivity was  the  study  of  methods  for  dealing  with  these  problems  in  the  context  of  BIE 
analysis.  The  results  of  this  study  are  discussed  in  Section  3. 
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SECTION  2 

ELASTIC  FRACTURE  MECHANICS  MODELING 
2.1  INTRODUCTION 

Gas  turbine  engine  disk  structures  are  often  life  limited  due  to  the  growth  of  fatigue  cracks. 
Such  cracks  have  been  found  to  fall  into  two  classes;  subsurface  or  buried  cracks  originating 
from  intrinsic  defects  and  surface  cracks  initiated  by  fatigue  loading  of  initially  defect-free 
structural  notches.  The  techniques  of  stress  analysis  and  fracture  mechanics  have  allowed  re- 
liable and  conservative  fatigue  life  prediction  for  buried  cracks,  and,  combined  with  appro- 
priate process  control  and  inspection  methods,  have  allowed  satisfactory  low  cycle  fatigue 
(LCF)  lives  to  be  achieved  for  structures  with  buried  defects.  The  LCF  life  prediction  prob- 
lem for  surface  cracks  is  considerably  more  complex.  Until  recently,  the  complexity  of 
both  the  crack  geometry  and  the  stress  field  near  structural  details,  such  as  rim  slots  and 
bolt  holes,  has  precluded  the  use  of  linear  elastic  fracture  mechanics  to  predict  propagation 
life  for  surface  cracks.  As  a result,  the  fatigue  life  of  structural  details  with  stress  concen- 
trations was  conservatively  estimated  by  predicting  the  number  of  cycles  to  initiate  a surface 
flaw  of  some  specified  size  and  ignoring  the  remaining  propagation  life. 


The  practical  extension  of  elastic  fracture  mechanics  techniques  to  the  prediction  of  propa- 
gation of  surface  (including  corner)  cracks  is  based  on  the  use  of  the  weight  function  tech- 
nique originally  developed  for  two-dimensional  problems  (Reference  1 7).  Development  of 
an  appropriate  geometrical  model  for  the  surface  crack  has  allowed  extension  of  the  weight 
function  technique  to  three-dimensional  problems  (Reference  18).  The  crack  modeling  stud- 
ies and  the  basic  improvements  in  three-dimensional  crack  modeling  capability  carried  out 
under  the  present  contract  have  allowed  the  development  at  Pratt  & Whitney  Aircraft  of  an 
extensive  stress  intensity  factor  data  base,  calibrated  with  surface  crack  growth  data. 

Further,  the  generally  improved  capability  for  calculation  of  stress  concentrations  in  com- 
plex structural  details  using  the  boundary-integral  equation  (BIE)  method,  which  has  been 
developed  both  under  this  contract  and  under  Pratt  & Whitney  Aircraft  in-house  programs, 
has  a direct  impact  on  both  the  initiation  and  propagation  problems  for  surface  flaws.  The 
system  for  the  prediction  of  LCF  cycles  to  surface  flaw  initiation  is  based  in  part  on  know- 
ledge of  the  local  concentrated  stress.  In  addition,  this  local  stress  field  in  the  uncracked 
part  is  used  in  the  weight  function  method  for  the  calculation  of  propagation  life. 

2.2  PREDICTION  OF  PROPAGATION  LIFE  FOR  SURFACE  FLAWS 

2.2.1  The  Weight  Function  Method 

It  has  been  shown  in  Reference  1 7 that  for  a linear-elastic  body  in  plane  strain  loaded  sym- 
metrically about  a crack,  the  Mode  1 stress-intensity  factor  for  any  load  system  can  be  cal- 
culated if  the  stress  distribution  in  the  uncracked  body  is  known  and  the  stress-intensity 
factor  and  displacement  field  are  known  for  a single  load  state.  In  particular,  in  the  absence 
of  body  forces 
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where  K’"*  and  t'"'  are  the  stress-intensity  factors  and  surface  tractions  of  Load  State  2 and 
u"  ’ is  the  displacement  field  for  Load  State  1 . H denotes  an  appropriate  elastic  modulus 
for  plane  stress  or  plane  strain,  and  C is  the  crack  length. 


In  the  case  mentioned  previously,  the  crack  is  one  dimensional,  and  its  stress  singularity  is 
defined  by  a single  stress-intensity  factor.  In  the  case  of  a crack  in  a three-dimensional 
geometry,  the  stress-intensity  factor  is  a function  of  position  on  the  crack  face.  In  princi- 
ple, the  stress-intensity  factor  for  any  load  state  can  be  determined  in  a manner  analogous 
to  Lq.  (b).  In  practice,  the  determination  of  the  required  information  on  variations  of  u'  * 
with  crack  geometry  is  too  difficult. 


In  order  to  exploit  the  weight-function  method  in  three-dimensional  problems,  it  is  nec- 
essary to  reduce  the  number  of  degrees  of  freedom  defining  the  crack  shape.  Service  and  ex- 
perimental experience  (References  6,  7,  and  19)  have  shown  tiiat  surface  and  corner  cracks 
tend  to  have  a part  elliptical  shape  of  rather  moderate  aspect  ratio  throughout  their  propa- 
gation life.  This  makes  plausible  an  approach  in  which  two  degrees  ot  Ireedom  (the  semiaxes 
a and  b)  are  allowed  to  determine  the  ellipitcal  crack  shape.  This  approach_( followed  in 
Reference  18)  allows  the  definition  of  the  averaged  stress  intensity  factors  and  Kj^,_which 
are  weighted  averages  of  the  stress  intensity  factor  distribution  along  the  crack  front.  and 
Kjj  are  then  related  by  the  material  crack  growth  model  to  independent  growth  of  the  two 
axes  of  the  elliptical  crack.  In  the  two  degree  of  freedom  model  this  procedure  leads  to  a 
coupled  pair  of  ordinary  differential  equations  for  crack  growth,  since  each  R depends  on 
the  current  shape  of  the  crack. 


Practical  application  of  this  two  degree  of  freedom  model  requires  an  efficient  method  for 
evaluating  and  R^  for  various  load  states  and  ellipse  aspect  ratios,  liquations  similar  to 
Hq.  6 allow  evaluation  of  and  Rj,  for  an  arbitrary  load  state  if  SW/QAl^ 
and  9W/9Al^  are  known  for  a single  reference  load  state.  W is  crack  opening  displacement 
and  A is  crack  area.  Calculation  of  these  quantities  requires  solution  of  three-dimensional 
fracture  mechanics  problems  and  is  too  expensive  for  routine  use.  As  an  alternative,  a data 
base  has  been  set  up,  using  three  dimensional  BIE  analysis,  in  which  are  stored  tlie  required 
data  for  buried,  surface  and  corner  cracks  in  a bar  of  rectangular  cross  section  for  a variety 
of  ellipse  aspect  ratios  and  bar  thicknesses.  Interpolation  using  this  data  allows  elficient 
calculation  of  and  for  arbitrary  aspect  ratios  and  thicknesses. 

2.2.2  Fracture  Mechanics  Calculations  for  Surface  Cracks 


Several  different  strategies  have  been  explored  in  recent  years  for  BlI:  modeling  of  cracked 
three-dimensional  structures.  The  more  general  approach  is  to  model  the  exterior  ot  the 
cracked  body  and  the  crack  surface,  but  not  the  unbroken  region  ahead  of  the  crack,  in 
which  the  traction  distribution  is  singular.  This  mo  leling  technique  is  applicable  without 
symmetry  restrictions  on  either  geometry  or  loading.  The  first  approach  along  this  line  was 
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the  idealization  of  the  crack  as  an  open  notch.  Useful  results  were  obtained  with  this  tech- 
nique (Reference  20),  but  it  poses  a major  dilemma.  The  accuracy  of  the  results  suffers  if 
the  notch  is  too  thick,  that  is  if  the  (coincident)  upper  and  lower  surfaces  of  the  mathemat- 
ical crack  are  modeled  too  far  apart.  On  the  other  hand,  the  equation  system  becomes 
badly  conditioned  if  the  surfaces  are  too  close  to  one  another. 

One  way  to  avoid  these  modeling  difficulties  is  the  development  of  a new  mathematical 
formulation  which  allows  modeling  of  both  crack  surfaces  in  the  same  plane  without  produc- 
ing a singular  BIE  coefficient  matrix.  During  the  first  phase  of  this  contract,  a new  BIE  was 
developed  for  three-dimensional  bodies  containing  a buried  plane  crack  of  arbitrary  shape 
(Reference  4).  The  displacement  interpolation  functions  of  the  only  three-dimensional  BIE 
code  (BINTEQ)  available  at  that  time  were  not  smooth  enough  to  allow  implementation  of 
the  new  BIE.  The  new  BIE,  in  conjunction  with  the  higher  order  BIE  codes  now  available 
(References  12  and  14)  can  provide  the  closest  analogy  in  three-dimensions  to  the  exact 
two-dimensional  mathematical  crack  modeling  capability  of  the  BIE/CRX  code  (Reference 
21). 

An  alternative  approach  to  three-dimensional  fracture  mechanics  modeling  proceeds  from 
the  observation  that,  in  gas  turbine  engine  structures,  fatigue  cracks  nonnally  grow  in  mode 
I,  normal  to  an  applied  tension  field.  Tliis  observation  suggests  using  the  crack  plane  as  a 
modeled  symmetry  plane,  avoiding  both  problems  associated  with  the  open  notch  model. 

This  modeling  strategy  was  studied  extensively  during  the  first  phase  of  the  contract.  Cal- 
culations were  carried  out,  using  the  BINTEQ  code,  for  buried,  surface,  and  corner  cracks 
modeled  as  ellipses  or  part -ellipses.  The  modeling  strategy  employed,  data  reduction  pro- 
cedures developed,  and  results  obtained  are  discussed  fully  in  References  4 and  7.  The  re- 
sults obtained  were  also  calibrated  against  experimentally  derived  mode  I stress-intensity 
factor  data  for  a growing  surface  flaw  (Reference  6).  Finally,  fracture  mechanics  data  ob- 
tained using  this  modeling  approach  has  been  used,  in  combination  with  the  weight  func- 
tion method,  to  calibrate  externally  generated  fatigue  life  data  for  corner  cracks  (References 
6 and  19). 

The  work  cited  above  has  been  shown  to  lead  to  satisfactory  predictions  of  fatigue  crack 
growth  for  design  system  purposes.  The  remaining  questions  center  in  two  areas; 

1.  Interaction  of  part  geometry  with  crack  growth.  The  improved  modeling  efficien- 
cy of  higher  order  BIE  codes  such  as  BASQUE  will  allow  resolution  of  such  effects. 

2.  Improved  efficiency  and  absolute  accuracy  of  BIE  calculations  in  the  near  crack  tip 
region.  This  characteristic  becomes  of  particular  importance  in  the  study  of  the 
crack/free  surface  intersection.  Work  conducted  as  part  of  the  current  contract 
has  been  devoted  to  the  development  and  evaluation  of  modified  BIE  crack  tip 
elements.  This  work  is  discussed  in  Section  2.3. 


J 


8 


PRATT  & WHITNEY  AIRCRAFT  GROUP 


2.3  USE  OF  BIE  SINGULARITY  EMEMENTS 
2.3. 1 The  Motlified  BIE  Crack  Tip  Element 


The  use  of  the  symmetrical  crack  modeling  strategy  for  mode  I loading  requires  the  approxi- 
mation of  near  crack  tip  displacements  and  tractions  in  terms  of  the  BIE  interpolation  func- 
tions. Considering  a straight  crack  front,  the  variation  of  both  displacement  and  traction  is 
linear  in  r (distance  from  the  crack  front)  in  BINTEQ.  In  the  BASQUE  code,  if  midpoint 
nodes  are  placed  at  the  geometrical  midpoints  of  the  element  sides,  the  boundary  data  varia- 
tion is  quadratic  in  r.  Since  the  near  crack  tip  variation  of  the  crack  opening  displacements 
is>/T  and  the  normal  traction  variation  ahead  of  the  crack  is  lA/T  , considerable  model 
refinement  is  necessary  to  achieve  acceptable  displacement  accuracy.  Since  singularities  are 
completely  absent  in  the  interpolation  functions,  the  tractions  near  the  crack  tip  can  never 
correctly  model  physical  behavior,  even  with  mesh  refinement. 

It  has  been  observed  that  proper  placement  of  certain  midpoint  nodes  in  three-dimensional 
isoparametric  finite  elements  leads  to  special  elements  with  the  appropriate  near  crack  tip 
displacement  and  traction  variation  (Reference  22).  A similar  element  has  been  developed 
for  the  higher  order  BIE  program.  The  midpoint  nodes  of  the  element  sides  normal  to  the 
crack  tip  are  repositioned  at  the  quarter  points  of  the  element  side,  relative  to  the  crack  tip. 
The  relationship  between  the  intrinsic  coordinate,  , and  r,  the  physical  distance  from  the 
crack  front,  is  then  (where  C = length  of  element  side) 


e 4 


(7) 


The  displacements  and  tractions  are  quadratic  in  J j . Using  the  shape  function  definitions  of 
Figure  2,  it  is  found  that  the  variation  of  Uj,  t^  in  the  physical  variable  is: 

( r— 

I t (r)|  " <8) 

In  particular,  the  crack  opening  displacement  Ujfr)  is 

u,(r)  = Aj  + Ajr  (9) 

since  Uj  (0)  = 0.  The  modified  BIE  crack  tip  element,  used  on  both  sides  of  the  crack,  thus 
allows  appropriate  variation  of  all  displacement  components  throughout  the  crack  tip  region. 
Since  boundary  displacements  and  tractions  are  independently  approximated  in  the  BIE 
method,  this  element  does  not  produce  the  \ k/T  traction  singularity  ahead  of  the  crack. 
This  is  in  contrast  to  the  finite  element  method  in  which  the  quarter  point  nodal  location 
automatically  gives  stresses  with  a 1/v^ singularity,  since  they  are  derived  by  differentiation 
of  the  displacement  interpolation  functions. 

file  modified  crack  tip  clement  has  been  used  for  two  studies;  a part-circular  surface  crack 
and  a center  cracked  test  specimen. 
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Figure  3 sliows  the  BASQUF  map  tor  the  eireular  eraek  problem.  One-eighth  of  the  body 
was  modeled.  The  v-^  plane  was  treated  as  an  unmodeled  symmetry  plane:  the  .\-z  plane 
was  modeled  to  allow  solution  of  the  buried  or  surface  crack  problem  by  a change  in  boundary 
conditions.  Figure  4 compares  the  crack  opening  displacements  for  the  buried  crack  for 
the  standard  and  modified  crack  tip  elements  using  both  linear  and  quadratic  boundary  data 
\ariation.  The  use  of  the  modified  crack  tip  element  is  seen  to  give  significant  improvement  | 

in  absolute  accuracy  of  displacements,  and  thus  K,  values,  for  the  same  map.  The  use  of  the 
modified  element  causes  no  increase  in  computer  run  time.  Figure  5 shows  the  variation  of 
K|  with  crack  front  location  for  the  surface  crack  which,  as  e.xpected,  confirms  results  re- 
ported earlier.  The  impact  of  the  crack  tip  element  is  not  a change  in  K,  variation  for  the 
surface  crack.  Rather,  it  is  that  the  improved  absolute  accuracy  allows  direct  K,  calculation 
without  resorting  to  the  rather  elaborate  scaling  methods  used  to  connect  BINTEQ  results 
to  a known  analytical  solution.  y 


Discussion  of  the  analysis  of  the  center  cracked  test  specimen  is  deferred  until  after  the  de- 
scription of  the  traction  singularity  element. 

2.3.2  The  Traction  Singularity  Element 


As  noted  above,  the  modified  crack  tip  element  does  not  allow  singular  tractions  at  the  crack 
tip.  It  was  found  possible  to  modify  the  BASQUE  code  to  produce  the  singularity  in  normal 
traction  on  the  first  element  ahead  of  the  crack.  The  shape  functions  for  traction  on  this 
element  were  modified  by  multiplication  by: 


0 tr) 


--  (1-v^) 
C 


(10) 


Because  of  programming  constraints,  the  modification  was  actually  made  in  the  routines 
which  evaluate  the  point  load  solutions.  The  function  0 (r)  is  O ( 1 A/  r)  for  small  values  of 
r and  is  one  at  r =C:  thus  the  singularity  is  produced  and  traction  continuity  is  guaranteed  on 
the  element  boundary  at  r = 8. 

Use  of  Eq.  (10)  with  a standard  element  gives 


(11) 


iM 


+ B,  v'T  + B,  r' 


while  its  use  with  the  modified  crack  tip  element  gives 

|t}  = — + B,  +B,V^ 

V ■■ 


(12) 


The  higher  order  terms  in  the  traction  expression.  Eq.  ( 1 1 ).  are  the  appropriate  ones. 
Unfortunately,  at  the  present  time  one  must  either  use  a crack  tip  element  ahead  of  the 
crack  and  accept  Fq.  ( I 2 1 or  use  a standard  element  and  obtain  Eq.  (II).  The  former  course 
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has  been  found  to  give  better  results  and  was  used  in  the  study  of  the  center  cracked  test 
specimen.  The  work  carried  out  to  date  has  shown  that  the  traction  singularity  element  has 
utility  which  justifies  the  programming  effort  required  to  remove  this  restriction  cited  above. 

2.3.3  The  Center  Cracked  Test  Specimen 

The  second  problem  analyzed  was  that  of  the  center  cracked  test  specimen  shown  in  Figure 

6.  One-eighth  of  the  geometry  was  modeled,  using  two  unmodeled  symmetry  planes  but 
modeling  the  crack  plane.  The  BIE  (BASQUE)  maps  used  in  the  analysis  are  shown  in  Figure 

7.  Both  two-dimensional  (plane  strain)  and  three-dimensional  analyses  were  carred  out  by 
changing  boundary  conditions  on  the  plane  y = t/2.  Singularity  finite  element  results  are 
available  for  the  three-dimensional  problem  (Reference  23).  Attention  was  given  to  several 
different  issues,  in  particular: 

1 . The  effect  of  the  use  of  crack  tip  and  traction  singularity  elements; 

2.  Evaluation  of  K,  without  extrapolation;  and 

3.  Exploration  of  the  region  near  the  crack  tip/free  surface  intersection  in  the  three- 
dimensional  problem. 

The  specimen  was  first  analyzed  in  plane  strain  using  Map  1 and  a variety  of  choices  of 
element  type  and  boundary  data  variation.  Figure  8 gives  the  stress  component  in  the 
direction  of  the  applied  load  for  the  various  cases  considered.  Interior  stresses  (0  = 90°) 
were  considered  as  it  was  found  that  they  were  more  accurate  than  the  traction  values  in 
the  crack  plane.  This  result  is  attributed  to  the  averaging  effect  of  the  boundary  integral 
representation  for  interior  stresses.  Further,  it  was  intended  to  use  stress,  rather  than  dis- 
placement, as  the  basis  for  K,  calculations,  avoiding  any  plane  stress  or  plane  strain  assump- 
tions in  the  the  three-dimensional  problem. 

Several  comments  can  be  made  about  the  results  of  Figure  8: 

1 . The  use  of  quadratic  variation  and  crack  tip  elements  gave  the  best  results; 

2.  The  divergence  of  the  finite  geometry  plane  strain  (BIE/CRX)  solution  from  the 
1 /v/~rapproximation  sets  an  outer  limit  (in  terms  of  r/a)  to  the  data  to  be  used 
in  calculating  K, ; and 

3.  The  boundary  layer  effect  of  finite  surface  element  size  on  the  three-dimensional 
BIE  interior  stress  calculation  sets  an  inner  limit  to  r/a. 

Based  on  these  observations,  it  was  decided  to  evaluate  K,  using  the  relationship 

o^(r,  0)  (13) 

with  0 = 90°,  0.01  < r/a  < 0.1 . Either  point  estimates  of  Kj  or  an  average  value  of  K, 
over  the  interval  (in  r/a)  can  be  used.  In  the  present  study  point  estimates  at  r/a  = 0.025 
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and  0.05  were  taken  because  a major  objective  in  the  three-dimensional  problem  was  to 
examine  stress  behavior  as  a function  of  r/a  near  the  free  surface. 

It  must  be  emphasized  that  K,  values  were  not  extrapolated  back  to  r/a  = 0.  Such  extrapola- 
tion was  felt  to  be  inappropriate,  at  least  in  the  BIE  analysis,  because  of  the  risk  of  misinter- 
preting the  numerical  boundary  layer  as  a real  variation  of  K,  with  r/a.  In  addition,  an  extra- 
polation technique  would  be  especially  dangerous  near  the  crack  tip/free  surface  intersection 
where  the  nature  of  the  stress  singularity  is  not  known.  Even  ignoring  numerical  boundary 
layer  effects,  the  basic  information  needed  to  guide  meaningful  extrapolation  does  not  yet 
exist. 

After  the  preliminary  work  done  with  the  relatively  crude  Map  1 , the  specimen  was  analyzed 
in  plane  strain  using  Map  2.  Modified  crack  tip  elements  were  used,  both  with  and  without 
the  traction  singularity.  Table  I contains  the  crack  opening  displacements  for  both  BASQUE 
analyses  at  three  different  locations  through  the  specimen  thickness.  The  results  from  a con- 
verged BIE/CRX  analysis  are  included  as  well.  There  is  significant  improvement  in  results 
near  the  crack  tip  due  to  inclusion  of  the  traction  singularity.  Further,  the  results  are  uni- 
form through  the  specimen  thickness,  indicating  that  the  BIE  map  on  the  plane  y = t/2  is 
adequately  refined. 


TABLE  1 

Ew/a  a : CRACK  OPENING  DISPLACEMENTS 
(PLANE  STRAIN) 


3D-Map  2 3D-Map  2 

crack  tip  element  crack  tip  element 

Plane  strain  with  traction  singularity  without  traction  singularity 


r/a 

(BIE/CRX) 

O 

II 

y/t  = 0.25 

y/t  = 0.5 

o 

II 

y/t-0.25 

If 

P 

0.0125 

0.399 

0.398 

- 

0.398 

0.381 

- 

0.380 

0.05 

0.795 

0.791 

0.792 

0.790 

0.784 

0.785 

0.785 

0.125 

1.243 

1.238 

- 

1.236 

1.234 

- 

1.233 

0.2 

1.552 

1.547 

1.547 

1.545 

1.543 

1.544 

1.541 

0.6 

2.444 

2.435 

- 

2.433 

2.432 

- 

2.429 

1.0 

2.693 

2.688 

2.688 

2.685 

2.684 

2.684 

2.682 

i 
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The  improvement  in  accuracy  due  to  the  traction  singularity  element  carries  over  to  the  K, 
calculation  as  well.  Table  II  compares  the  K,  values  obtained  with  and  without  the  singularity 
element  at  r/a  = 0.05.  The  BIE/CRX  result  given  is  believed  accurate  to  within  0.25  percent. 

TABLE  II 

Kj/<V^  FOR  CENTER  CRACKED  PANEL' 


(PLANE  STRAIN) 

Analysis 

From  Interior 

Stress  (9  = 90°) 

From  Crack 

Opening  Displacements 

BIE/CRX 

1.42)2 

1.421 2 

BASQUE  — with 
traction  singularity 

1.420 

1.416 

BASQUE  - without 
traction  singularity 

1.436 

1.395 

1 - evaluated  at  r/a  = 0.05 

2 - evaluated  using  path  independent  integral.  Reference  21. 


It  is  clear  that  the  use  of  the  traction  singularity  element  improves  both  the  accuracy  and 
the  consistency  between  the  stress  and  displacement  fields. 

Solution  of  the  three-dimensional  problem  was  carried  out  using  Map  2;  later,  the  more 
refined  Map  3 was  also  used  with  no  significant  change  in  results.  Figure  9 gives  nondimen- 
sional  crack  opening  displacements  for  the  three-dimensional  problem  at  both  the  center 
and  free  surface  of  the  specimen.  Results  are  given  for  the  BASQUE  analysis,  the  singular 
finite  element  analysis  of  Raju  and  Newman  (Reference  23),  and  for  comparison,  for  the 
BIE/CRX  solution.  The  results  of  the  BASQUE  analysis  show  displacements  significantly 
larger  than  the  finite  element  results  at  both  the  center  and  free  surface  of  the  specimen. 

K,  at  the  center  of  the  specimen  is  within  1 percent  of  the  plane  strain  value. 

The  effect  of  the  free  surface  is  shown  in  Figure  1 0.  Stress  is  plotted  as  a function  of  posi- 
tion through  the  specimen  thickness.  The  BASQUE  results  are  given  at  r/a  = 0.025  and  0.05 
The  finite  element  results  are  based  on  extrapolation.  The  BIE  and  finite  element  results 
agree  qualitatively,  showing  an  increase  in  stress,  peaking  at  y/t  ^ 0.42,  followed  by  a rapid 
drop  off  in  a boundary  layer  near  the  free  surface.  The  BASQUE  stress  data  (at  both  r/a 
locations)  rise  to  a higher  peak  {+  2.5  percent)  closer  to  the  free  surface  than  the  finite  ele- 
ment data.  The  BASQUE  stresses  then  drop  faster  to  lower  values  at  the  free  surface,  rela- 
tive to  the  singular  finite  element  data. 


In  addition,  comparison  of  the  BIE  results  at  the  two  different  r/a  values  indicates  a possible 
change  in  the  nature  of  the  stress  singularity  at  the  free  surface.  The  behavior  of  stress  near 
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the  peak  stress  location  is  consistent  with  a singularity  of  higher  order  than  1 'y/T"  It  is 
reasonable  to  suppose,  however,  that  beliavior  near  the  crack  tip/free  surface  intersection 
depends  on  position  relative  to  botli  the  crack  and  the  free  surface.  It  may  well  be  inappro- 
priate to  interpret  the  data  near  the  free  surface  simply  in  terms  of  distance  from  the  crack 
tip. 

Finally,  it  is  possible  that  the  nature  of  behavior  near  a free  surface/crack  tip  intersection  can- 
not be  resolved  by  numerical  modeling.  In  the  work  discussed  here,  it  has  been  verified  that 
the  behavior  shown  was  not  induced  by  the  use  of  the  traction  singularity  element.  Never- 
theless. both  the  BIF  and  finite  element  methods  assume  specific  interpolation  functions, 
and  the.se  assumptions  limit  the  solution  behavior  which  can  be  e.xhibited  without  excessive 
mesh  refinement.  Pending  the  analytical  resolution  of  the  nature  of  the  singularity  (if  any) 
at  the  crack  tip'  tree  surface  intersection,  the  most  promising  approach  is  the  implementa- 
tion of  the  Oat  crack  BIE.  This  approach  would  completely  eliminate  the  traction  modeling 
problem  and  allow  needed  mesh  refinement  to  be  concentrated  on  the  more  accurate  defini- 
tion of  displacements  on  the  free  surface  and  crack  near  their  intersection. 


( 
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SECTION  3 

BIE  ANALYSIS  FOR  INHOMOGENEOUS  MATERIALS 

3.1  INTRODUCTION 

In  gas  turbine  engine  applications,  inhomogeneous  elastic  material  properties  are  usually 
encountered  as  a result  of  spatial  temperature  variation  and  temperature  dependence  of  elas- 
tic material  properties.  This  situation  is  particularly  acute  in  the  combustor  and  turbine  sec- 
tions of  an  engine.  In  these  sections,  BIE  analysis  can  profitably  be  applied  to  turbine  blade 
attachments  and  the  turbine  disk  rim  region  using  a three-dimensional  code  and  to  axisym- 
metric  analysis  of  a disk  using  a BIE  code  such  as  that  described  in  Reference  24. 

The  analysis  required  is  a mission  analysis;  that  is,  stress  analysis  of  a part  must  be  carried 
out  at  several  points  of  the  engine  operating  cycle.  The  temperature  distribution  in  the  part 
is  neither  uniform  nor  steady  state  at  most  operating  points;  thus,  effective  use  of  the  BIE 
method  requires  development  of  both  strategy  for  efficiently  handling  the  transient  thermo- 
elastic problem  and  a technique  for  treatment  of  material  inhomogeneity.  Tlie  remainder  of 
this  section  discusses  possible  techniques  for  the  introduction  of  inhomogeneous  elastic 
material  properties  in  BIE  analyses  and  closes  by  proposing  a strategy  which  fits  the  require- 
ments of  a mission  analysis. 

3.2  APPROXIMATE  BIE  ANALYSIS  FOR  INHOMOGENEOUS  MATERIALS 

Even  if  elastic  material  properties  are  assumed  to  vary  in  space,  Eq.  ( 1 ) and  (2)  can  still  be 
formally  derived.  It  is  not,  however,  possible  to  derive  an  analytical  expression  for  the  point 
load  functions  U^  (p,  Q),  Tj^tp,  Q)  for  a general  inhomogeneous  material,  even  if  it  is  assumed 
to  be  isotropic.  Tlie  absence  of  such  an  expression  precludes  the  evaluation  of  the  point  load 
functions  by  a numerical  technique  such  as  that  described  in  the  Appendix  for  homogen- 
eous anisotropic  materials.  There  remain  three  approaches:  derivation  of  point  load  .solu- 
tions for  special  types  of  inhomogeneity,  material  property  averaging,  and  iterative  techniques. 
These  approaches  are  explored  in  the  paragraphs  below. 

The  first  possibility,  the  construction  of  point  load  functions  for  special  type.s  of  mhomogeneify 
is  attractive  in  some  situations  but  not  in  the  gas  turbine  engine  environment  It  is  to  be  expected 
that  such  special  point  load  solutions  could  he  derived  only  for  very  simple  material  pro- 
perty variations,  such  as  linear  variation  in  the  spatial  variables.  This  means  that,  in  the  so- 
lution of  any  particular  problem,  the  actual  material  property  variation  must  be  approxim- 
ated using  those  variations  for  which  special  point  load  functions  exist.  In  gas  turbine  en- 
gine parts  the  material  property  variation  is  not  basically  spatial,  rather  it  is  due  to  temper- 
ature dependence  of  the  elastic  material  properties  and  the  existance  of  a highly  nonuniform 
temperature  field.  This  environment  leads  to  a complicated  spatial  variation  of  material  pro- 
perties and  makes  the  required  approximation  difficult  if  not  impossible.  This  approach  is 
therefore  not  appropriate  for  the  structural  analysis  of  gas  turbine  engine  parts  althougli  it 
could,  at  some  time,  prove  useful  in  other  applications. 

A second  approach  is  the  approximation  of  the  inhomogeneous  material  by  a homogeneous 
material,  using  average  elastic  material  properties  for  the  part.  Tlie  calculation  of  the  average 
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properties  poses  no  problem,  as  the  normal  design  process  requires  thermal  analysis  of  the 
part  prior  to  stress  analysis.  The  results  of  an  in-house  study  carried  out  using  an  axisym- 
metric  finite  element  code  are  discussed  below,  and  these  results  clearly  demonstrate  the 
feasibility  of  the  approach. 

The  problem  studied  was  the  stress  analysis  of  a typical  turbine  disk  under  combined  cen- 
trifugal, thermal,  and  blade  loads.  The  cross-section  of  the  disk  is  shown  in  Figure  1 1 along 
with  isotherms  from  the  thermal  analysis  of  the  disk.  The  analysis  was  carried  out  under  two 
different  material  property  assumptions; 

1.  E and  a (Young’s  modulus  and  coefficient  of  thermal  expansion)  evaluated  based 
on  temperature  within  each  element  (Table  HI). 

2.  E evaluated  at  the  average  disk  temperature  (567°F)  and  a at  its  local  value  in 
each  element. 

Density  and  Poisson’s  ratio  were  held  constant,  although  that  is  not  a necessary  restriction. 
The  results  of  the  two  analyses  are  compared  for  the  disk  bore  and  disk  rim  in  Tables  IV 
and  V,  respectively.  The  locations  at  which  the  comparisons  were  made  are  shown  in  Figure 
II. 

TABLE  111 


TEMPERATURE  DEPENDENCE  OF  MATERIAL  PROPERTIES 


Temperature 

E 

a 

(°F) 

(psi) 

(in/in/°F) 

0 

31.25x10^ 

6.76x10'^ 

100 

30.75 

6.895 

200 

30.30 

7.01 

300 

29.90 

7.13 

400 

29.50 

7.25 

500 

29.09 

7.36 

600 

28.62 

7.46 

700 

28.19 

7.57 

800 

27.70 

7.67 

900 

27.20 

7.78 

1000 

26.70 

7.88 

1100 

26.19 

7.99 

1200 

25.60 

8.11 

1300 

25.00 

8.25 

1400 

24.40 

8.40 

1500 

23.70 

8.62 

J 


PRATT  & WHITNEY  AIRCRAFT  GROUP 


TABLE  IV 

EFFECT  OF  MATERIAL  PROPERTY  AVERAGING  AT  DISK  BORE 


Location 


Ba 

B3 

B4 

B5 

Be 


Radial  Displacement  (in) 

Hoop  Stress  (psi) 

Analysis  1 

Analysis  2 

Analysis  1 

Analysis  2 

0.01522 

0.01526 

56,694 

58,465 

0.01675 

0.01679 

82,615 

84,604 

0.01833 

0.01838 

101,069 

103,142 

0.01922 

0.01925 

103,454 

105,840 

0.01803 

0.01804 

84,735 

86,663 

0.01647 

0.01647 

49,759 

51,560 

TABLE  V 

EFFECT  OF  MATERIAL  PROPERTY  AVERAGING  AT  DISK  RIM 


Radial  Displacement  (in.) 

Radial  Stress  (psi) 

Hoop  Stress  (psi) 

Location 

Analysis  1 

Analysis  2 

Analysis  1 

Analysis  2 

Analysis  1 

Analysis  2 

0.05583 

0.05570 

50,424 

50,453 



B2 

0.05624 

0.05613 

50,639 

50,644 

— 

— 

R3 

0.05639 

0.05629 

51,533 

51,596 

- 

- 

R4 

0.05636 

0.05627 

51,545 

51,586 



_ 

B5 

0.05617 

0.05607 

50,629 

50,618 

— 

— 

Be 

0.05573 

0.05563 

50,169 

50,212 

- 

- 

R7 







-17,337 

-18,601 

Bs 

- 

- 

— 

— 

-1,818 

-2,039 

B9 

- 

— 

- 

- 

31,687 

32,578 

Bio 

_ 

— 

— 



26,297 

27,155 

B„ 

- 

- 

- 

— 

14,480 

15,023 

Bu 

- 

- 

- 

- 

-16,300 

-17,022 

The  good  agreement  between  the  two  analyses  demonstrates  the  potential  of  the  BIE  method 
in  dealing  with  thermally  induced  material  property  inhomogeneity.  It  should  be  noted  that 
local  evaluation  of  a is  acceptable  in  the  BIE  method  since  it  does  not  enter  into  the  point 
load  functions. 

Substructuring  could  also  be  used  to  enhance  the  accuracy  of  the  material  property  averaging 
technique.  Reference  to  Figure  1 1 shows  clearly  that  decomposition  of  the  disk  cross  section 
into  two  or  three  appropriately  chosen  subregions  can  substantially  reduce  the  difference  be- 
tween extreme  and  average  temperature,  and  should  yield  a corresponding  increase  in  accuracy 
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l inally.  an  iteratJvc  tecliniijue  for  the  inhomogeneous  problem  has  been  shown  to  converge  in 
two-dimensional  analysis  ( Ret'erence  25).  To  summarize  the  results,  consider  the  equilibrium 
equation: 


%’i  <>4) 

In  the  problems  considered  here,  the  will  normally  be  the  centrifugal  and  thermal  body 
forces;  Oy  is  the  stress  tensor.  In  the  inhomogeneous  problem,  the  termsOj|.|  involve  non- 
vanishing space  derivatives  of  both  strains  (as  in  the  homogeneous  case)  and  material  pro- 
perties. Substituting  in  Tq.  (14)  in  terms  of  displacements  results  in: 


1 


+ —X  -i-  F.  = 0 

^ I j 


(15) 


where  F^  depends  on  derivatives  of  Uj  and  on  the  material  properties  (/a,  i>).  For  a homogeneous 
material,  F.  = 0,  and  F’q.  ( 1 5)  is  Navier’s  equation.  The  basic  result  of  Reference  25,  proved 
by  converting  Fq.  ( 15)  to  an  integral  equation,  is  that  Eq.  1 5 can  be  solved  by  an  iterative 
technique  whose  lowest  order  solution  is  that  for  a homogeneous  material.  Thus,  u°  is  the 
solution  of  Fq.  (15)  with  F°=  0.  Then  F^'  is  calculated  using  u°  and  u.'  is  determined  as  the 
solution  of: 


■ 1 ' I 


1 1 
+ - X + F 
/a  ' j 


= 0 


(16) 


Building  on  this  result,  a technique  for  extending  the  BIE  method  to  inhomogeneous  ma- 
terials is  to  first  solve  the  problem  with  homogeneous  (average)  material  properties  and  then 
improve  this  result  by  calculating  the  term  F^  and  applying  it  as  an  additional  body  force 
term  in  a second  BIE  solution.  The  results  cited  for  the  averaging  technique  alone  make  it 
unlikely  that  more  than  one  or  two  iterations  would  be  required  to  achieve  acceptable  ac- 
curacy in  practical  problems. 

3.3  A COMPUTATIONAL  STRATEGY  FOR  TRANSIENT  THERMAL  ANALYSIS 

As  was  mentioned  in  Section  1 , the  presence  of  body  forces,  f(q),  requires  the  introduction 
of  a volume  integral  for  each  source  point  P (x), 

fjj  (P,q)  f(q)  dv  (17) 


in  the  BIE.  Eq.  (2).  The  centrifugal  load  terms  can  be  converted  directly  to  a surface  integral, 
but.  for  a time  dependent  temperature  field,  the  thermal  loading  term  cannot  be  converted. 
Since  time  dependent  temperature  fields  are  characteristic  of  mission  analysis,  an  efficient 
method  will  be  required  for  the  evaluation  of  Eq.  (17).  Preliminary  studies  indicate  that 
the  evaluation  can  be  accomplished  without  undue  impact  on  Bll-  solution  time. 
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Heat  transfer  analysis  for  gas  turbine  engine  parts  is  typically  conducted  by  using  methods 
which  subdivide  the  region  into  a finite  number  of  two-  or  three-  dimensional  elements,  B^, 
each  of  which  has  a single  assigned  value  of  temperatue,  T,^ . Heat  transfer  analysis  is  norm- 
ally carried  out  prior  to  stress  analysis  so  that  the  geometrical  subdivision  of  the  part  and  the 
time  dependent  temperatures  are  known  to  the  stress  analyst.  The  use  of  the  explicit  expres- 


sion for  the  point  load  function  U y allows  Eq.  (17)  to  be  rewritten  (see  Reference  4)  as; 


i 


F(n,  E,  V)  G(P,  q)  T (q)  dv 


(18) 


where  a is  the  coefficient  of  thermal  expansion,  E is  Young’s  modulus,  v is  Poisson’s  ratio, 
T (q)  is  the  temperature  field  and  G (P,q)  is  a known  function  involving  derivatives  of  the 
point  load  function  U...  The  expression  can  then  be  expanded  as: 


n)  G(P,q)  T^(q)  dv 


Since  temperature  is  constant  within  each  element,  Eq.  (19)  can  be  evaluated  as; 
M 


(19) 


1 


^’k)T^  / G(P,q)dv 


CO) 


The  calculation  of  Eq.  (20).  using  a numerical  volume  integration  technique,  is  expensive, but 
it  is  required  only  once  since  the  individual  integrals  do  not  involve  temperature.  Once  the  in- 
tegrals are  available,  they  can  he  used  repeatedly  to  generate  the  thermal  load  vector  for  as 
many  time  points  as  desired.  Further,  incorporation  of  temperature  dependence  in  a is  straight- 
forward. since  it  occurs  outside  the  integral  sign  in  Eq.  ( 20).  The  temperature  dependence  of  E 
can  be  handled  using  the  averaging  technique  described  in  Section  3.2. 

Preliminary  estimates  indicate  that  the  overhead  cost  for  including  transient  thermal  analysis 
should  be  50  to  100  percent  relative  to  the  same  problem  under  mechanical  loads  alone,  but 
that  the  added  cost  for  each  additional  time  point  should  be  under  5 percent.  This  cost  makes 
the  technique  attractive  relative  to  finite  element  methods. 
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Figure  2 Quadratic  Isoparametric  Shape  Functions 


Figure  3 Boundary-Integral  Fquation  Map  for  a Circular  Crack 


Figure  4 Displacement  Variation  for  a Buried  Circular  Crack 
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Figure  1 1 Radial  Section  of  Turbine  Disk 
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SUMMARY 


The  boundary-integral  equation  method  is  particularly  well  suited  for  solution  of  stress 
concentration  and  elastic  fracture  mechanics  problems.  The  method  was  not  previously 
applicable  to  anisotropic  three  dimensional  problems  because  no  efficient  technique  existed 
for  calculation  of  the  required  point  load  solution  for  an  infinite  body.  A technique  has 
been  developed  to  evaluate  numerically  the  anisotropic  point  load  solutions,  and  used  to 
generate  data  bases  for  various  materials.  An  interpolation  technique  is  used  to  evaluate 
the  point  load  solutions  efficiently  within  a higher  order  boundary-integral  equation  code. 

‘ The  effectiveness  of  the  technique  is  verified  by  solution  of  problems  involving  both  uni- 

! axial  stress  states  and  stress  concentrations. 
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INTRODUCTION 


Increasing  structural  use  is  being  made  of  materials  with  anisotropic  elastic  material  pro- 
perties. Composites  have  been  used  with  substantial  weight  advantage  in  various  airframe 
structures.  Eutectic  and  directionally  solidified  alloys  are  finding  increasing  use  in  advanced 
gas  turbine  engines  to  provide  increased  strength  without  weight  and  performance  penalties. 

The  elastic  stress  analysis  required  for  such  anisotropic  materials  falls  generally  into  one  of 
two  classes,  the  general  analysis  of  an  entire  structure  (for  example,  a disk  or  turbine  vane) 
or  the  detailed  analysis  of  a critical  subsection.  The  general  stress  analysis  of  a structure 
does  not  usually  require  a precise  definition  of  high  stress  gradients  and  can  be  effectively 
carried  out  (in  both  two  and  three  dimensions)  using  generally-available  stress  analysis  com- 
puter codes. 

A detailed  and  accurate  calculation  of  stresses  in  regions  of  rapid  stress  variation  is  often 
required,  both  for  general  design  purposes  and  as  input  for  fatigue  life  calculations.'  Three- 
dimensional  analysis  is  often  needed  to  account  properly  for  the  effects  of  part  geometry 
and  loading.  The  boundary-integral  equation  method  is  particularly  well  suited  to  problems 
requiring  the  resolution  of  high  stress  gradients.’  In  two  dimensions,  it  already  provides  a 
numerical  technique  applicable  to  either  isotropic  or  anisotropic  materials.^  This  paper  de- 
scribes the  extension  of  the  technique  to  anisotropic  three-dimensional  elasticity.  A method 

for  the  numerical  calculation  of  the  anisotropic  point  load  solutions  is  presented  and  results 
of  boundary-integral  equation  analyses  for  several  anisotropic  materials  are  discussed. 

ANALYTICAL  FORMULATION 

The  Boundary-Integral  Equation  Method 

This  paper  does  not  attempt  a detailed  review  of  the  boundary-integral  equation  method. 
For  this  the  reader  is  referred  to  Cruse*  or  Lachat  and  Watson."*  We  simply  recall  that  the 
method  is  based  on  knowledge  of  the  point  load  solution  Ujj(x,  y)  for  an  infinite  body, 
which  gives  the  displacements  at  the  field  point  y due  to  a point  force  applied  at  the  source 
point  x^.  Use  of  the  reciprocal  work  theorem  and  appropriate  limit  operations  give  the 
boundary-integral  equation 


r 


(1)  + ^Tjj(x,  Y)uj(y)ds(y)  = ^yUjj(x,  y)tj(y)ds(y) 


where  Tjj  is  the  traction  point  load  solution  derived  from  Ujj  and  Uj,  tj  are  the  boundary 
data  on  the  region  S.  Appropriate  numerical  modeling  of  Uj,  tj  and  S then  reduces  (1 ) to  a 
set  of  linear  algebraic  equations.  One  of  the  major  advantages  of  the  method,  apparent  from 
(1),  is  that  no  interior  idealization  of  the  body  is  required.  Stresses  and  displacements  at 
selected  interior  points  can  be  evaluated  by  surface  integration  after  the  boundary  solution 
is  completed. 

Anisotropic  Point  Load  Solutions 


For  three-dimensional  isotropic  elasticity  the  point  load  function  U|j(x,  y)  is  the  well  known 
Kelvin  solution. 


(2a)  Uj  (x,  y)  = 


) (3-4r^)  g _1 (Xj-yj)  (xj-yj) 

47rplxyl  I 4(1-1^)  'i  '''  4(1-1')  |x-y|  |x-y( 


For  a general  anisotropic  material  the  point  load  solution  can  be  represented®'  ® as 


(2b)  U,j(x,y)  = 


Srr'^lxyl 


Kjjl 


The  line  integral  is  taken  on  the  unit  circle  in  the  plane  normal  to  (x-y)  and  passing  through 
X.  The  function  is 


(3)  Kjj'(t)  = 


] ^ 

^ijkm^k^  m 


where  the  Cjj|^f^  are  the  elastic  constants  of  the  material.  Representations  for  the  traction 
point  load  solution  Tjj  can  be  derived  from  (2a)  or  (2b).  The  representation  (2b)  of  Ujj,  and  the 
associated  Tjj,  are  not  directly  computable.  The  extension  of  the  boundary-integral  equation 
method  as  an  effective  numerical  technique  for  three-dimensional  anisotropic  problems  re 
quires  an  accurate  and  efficient  means  for  calculating  the  anisotropic  point  load  solutions. 
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EVALUATION  OF  ANISOTROPIC  POINT  LOAD  SOLUTIONS 
Reduction  of  Point  Load  Solutions  to  Computable  Form 

The  representation  (2b)  cannot  generally  be  evaluated  in  closed  form.  For  an  isotropic 

material  it  reduces  analytically  to  the  well-known  Kelvin  form.^  The  only  other  case  for 
which  a closed-form  solution  is  available  is  that  of  a transversely  isotropic  material.®  ’ ' " 

A second  approach  involves  series  expansions  of  (2b)  and  its  associated  traction  solution.  A 
uniformly  convergent  expansion  has  been  obtained"  , but  is  not  suitable  for  extensive 
computation. 

The  remaining  alternative  is  a direct  numerical  evaluation  of  the  functions,  Ujj  and  T-. 

One  approach  to  this  problem  is  indicated  by  Vogel  and  Rizzo’,  but  the  method  discussed 
is  quite  complex,  especially  for  the  function  T-  and  would  be  too  time  consuming  for  rou 
tine  numerical  use. 

The  present  approach  to  the  calculation  of  Ujj  and  Tjj  is  based  on  defining  the  modulation 
function. 


(4) 


Gij(ei,U2) 


(?)ds 


1^1=  1 

where  v^,  V2  define  the  orientation  of  the  vector2<-y-  Then 


(5) 


^ij 


— O G 

STT^Ix-yl  'I 


As  was  observed  in’  all  the  singular  behavior  of  Ujj  occurs  in  the  first  factor;  further  Gjj  is 
independent  of  Ix-y) . The  evaluation  of  Tjj  (the  traction  point  load  solution)  using  Hooke's 
Law  requires  the  calculation  of  the  displacement  gradients  Ujj  |^.  They  can  be  expressed  as 


Ujj,  k 


G-.  + 

87r2|x-y|3  '*  Srr^lx.yl 


A-.S 


*^ij,  a'^a,  k 
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The  known  (x-yC^  singularity  occurs  in  easily  evaluated  closed  form  expressions,  including 
the  derivatives  k-  The  derivatives  Gjj,  ^ are  all  non-singular,  with  the  exception  of  possi 
ble  removable  singularities  induced  by  the  choice  of  the  u^,  V2  coordinate  system.  The 
higher  order  derivatives  needed  for  interior  stress  calculations  can  be  similarly  expressed. 


The  modulation  function  Gjj  can  be  calculated  by  a straightforward  numerical  evaluation 
of  the  line  integral{3).  The  integrand  can  be  expressed,  in  closed  form,  as  a rational  func 
tion  of  degree  -2.  The  key  to  the  evaluation  of  the  derivatives  Gjj,  ^ is  that  the  integral  in 
(3)  can  be  transformed  so  that  the  integration  path  is  independent  of  the  orientation  of 
x-y;  as  a result  the  integrand  will  then  depend  explicitly  on  u-j  and  V2-  The  required  deriva 
tives  can  then  be  calculated  as 


'51  = 1 


by  differentiation  through  the  integral  sign.  Transformation  of  the  material  constants 
Cjjk)  is  not  required,  rather  an  explicit  substitution  ^j  = fj  (r?i , ^2)  is  carried  out.  Further 

since 


Kii  K 1 

')  jk 


6 


ik 


implicit  differentiation  gives 

*^ij  ' Q “ *^ik  *^kl'  a **'lj 

and  the  only  explicit  differentiation  required  is  that  of  the  quadratic  forms  (in  rjj,  172)  which 
are  the  elements  of  Kjj;  Any  type  of  numerical  differentiation  is  completely  avoided. 

Numerical  Evaluation  of  the  Modulation  Function 

The  calculations  described  have  been  carried  out  for  a variety  of  materials.  The  line  inte- 
grals were  evaluated  using  Simpson's  rule.  Table  1 shows  the  convergence  of  the  modula 
tion  function  as  the  integration  is  refined. 


,\-b 


TABLE  1 


Convergence  of  Modulation  Function* 

Points  on  Contour  G-|i  G^3 


8 

3.462192  X 10'^ 

3.959000  X 10'^ 

16 

3.664519  X 10'^ 

2.933159  X 10'^ 

32 

3.598752  X 10'^ 

2.959600  X 10  ® 

64 

3.600230  X lO'^ 

2.958070  X 10  ® 

128 

3.600341  X 10'^ 

2.958069  X 10  ® 

*0  = I'l  = 22,5°,  0 = r’2  = 67.5°;  face  centered  cubic  material 

In  order  to  test  the  algorithm  for  calculating  Tjj,  the  tractions  for  a face  centered  cubic  ma- 
terial were  integrated  over  the  surface  of  a sphere  containing  the  source  point.  The  results, 
shown  in  Table  2,  demonstrate  the  recovery  of  the  applied  load. 

TABLE  2 


Resultant  Load  on  Spherical  Surface* 


Applied  Load 

Resultant  Load  Direction 

Direction 

X 

Y 

Z 

X 

.999897 

.000048 

.000099 

Y 

.000030 

.999808 

-.000147 

Z 

-.000004 

.000004 

.999706 

*Unit  sphere  centered  at  (.1,  .2,  .3);  load  point  at  (0,  0,  0).  Face  centered  cubic 
material.  Numerical  integrations:  32  points  for  line  integrals. 
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The  technique  outlined  above  was  incorporated  in  a higher  order  three-dimensional  boundary 
integral  equation  code.'*  A few  very  simple  problems  were  solved  to  verify  the  method,  but 
relative  computing  times  were  so  large  that  the  method  was  unsuitable  for  any  realistic  prob- 
lem. A further  development  of  the  numerical  technique  to  achieve  practical  computing  times 
is  described  below. 

Application  to  Large  Scale  Stress  Analysis 

The  analytical  development  and  numerical  results  discussed  above  show  clearly  that  the 
anisotropic  point  load  solutions  can  be  computed  with  essentially  arbitrary  accuracy.  The 
fact  that  a face  centered  cubic  material  was  chosen  for  the  initial  investigation  implies  no 
restriction.  The  present  effort  also  treats  transversely  isotropic  and  orthotropic  materials; 
further,  the  extension  to  full  anisotropy  requires  only  the  algebra  to  calculate  the  explicit 
forms  of  the  integrands  in  (3)  and  (6). 

The  crucial  issue  in  determining  the  applicability  of  the  boundary-integral  equation  method 
to  anisotropic  problems  of  engineering  interest  is  that  of  efficiency.  The  lack  of  closed  form 
expressions  for  the  anisotropic  kernel  functions  requires  that  the  kernel  function  - boundary 
data  integrals  in  (1)  be  evaluated  numerically.  These  integrations  are  typically  the  most  time 
consuming  part  of  a boundary-integral  equation  analysis.  Using  the  closed  form  isotropic 
kernels,  the  numerical  integrations  in  the  boundary-integral  equation  code  used  require  an 
average  of  .0015  cpu  sec.  for  each  integration  point,  of  which  only  .0001  cpu  sec.  is  required 
for  the  actual  kernel  function  calculation.  By  contrast,  a single  calculation  of  Uj.  and  T.^  for 
an  anisotropic  material  by  direct  integration  (using  64  integration  points)  requires  .025  cpu 
sec.  This  implies  an  increase  in  overall  computing  times  by  a factor  of  12  to  16,  and  is 
clearly  unacceptable  for  any  practical  problem. 

The  solution  to  this  problem  is  the  substitution  of  an  interpolation  technique  for  direct 
numerical  integration  in  the  evaluation  of  U,.  and  T^^.  The  technique  is  suggested  by  the 
fact  that  G,.  and  its  derivatives  are  all  smooth  functions  of  the  variables  v,,  v^.  The  method 
has  been  implemented  and  tested  for  a variety  of  materials. 

G||  and  its  derivatives  were  expressed  as  functions  of  </>  and  6 (the  two  angles  in  spherical 
coordinates)  in  order  that  interpolation  could  be  carried  out  in  a rectangular  table 
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(0  < 7r,  0 < 6 < 7r).  Table  entries  were  calculated  by  direct  numerical  integration  using 

64  integration  points  for  the  line  integral  evaluations.  Interpolation  in  the  tables  was  carried 
out  using  quadratic  or  cubic  Lagrange  interpolation. 

To  verify  the  interpolation  technique,  the  tractions  were  again  numerically  integrated  over 
a sphere.  The  recovery  of  the  applied  unit  loads  is  shown  in  Table  3.  For  these  tests  a 17  x 33 
table  was  used  for  G...  The  accuracy  obtained  using  the  interpolated  point  load  solution  is 
essentially  the  same  as  that  obtained  using  the  closed  form  isotropic  kernels.  In  addition, 
the  load  recovery  for  the  anisotropic  (transversely  isotropic)  material  is  equivalent  to  that 
for  an  isotropic  material.  The  overall  level  of  accuracy  is  slightly  less  than  shown  in  Table  2 
because  a coarser  integration  mesh  was  used  on  the  sphere. 

TABLE  3 

Resultant  Load  on  Spherical  Surface 


Transversely 


i 

j 

Applied 

Load 

Isotropic  - 
Exact  T|| 

Isotropic  - 
Quadratic 
Interpolation 

Isotropic  - 
Cubic 

Interpolation 

Isotropic 

Quadratic 

Interpolation 

1 

1 

1 

1.0 

.99860 

.99794 

.99845 

.99874 

2 

0 

.00007 

-.00007 

.00005 

-.00017 

3 

0 

-.00017 

-.00007 

-.00017 

-.00002 

2 

1 

0 

.00007 

-.00014 

-.00005 

-.00015 

2 

1.0 

.99870 

.99818 

.99837 

.99917 

3 

0 

-.00034 

-.00014 

-.00034 

-.00037 

3 

1 

0 

-.00020 

-.00010 

-.00020 

-.00009 

2 

0 

-.00040 

-.00020 

-.00040 

-.00019 

3 

1.0 

1.00069 

1.00150 

1.00101 

1.00102 

Modulation  function  data  bases  were  generated  for  a variety  of  materials  and  stored  on  mag- 
netic tape  for  use  by  a suitably  modified  version  of  the  boundary-integral  equation  program. 
The  table  size  used  was  33  x 65.  The  time  required  for  the  data  base  generation  was  15  cpu 
min.  for  each  material,  including  the  calculation  of  all  the  derivatives  Gjj,  which  would  be 
required  for  interior  stress  calculations.  The  material  constants  used  are  shown  in  Table  4. 


TABLE  4 


Material  Constants  for  Data  Base  Generation 


Constant* 

1 - Isotropic 

E = 18.07  X 106  psi 
1^=  .38931 

2 - Isotropic 

E = 18x  106  psi 

V = .3 

3-Transve'’sely  Isotropic 
CobaltB 

C11 

36  X 106 

24.23  x 106 

44.52  X 106 

C22 

= Cli 

= Cii 

= Cn 

C33 

= Cii 

= Cii 

51.94  X 106 

C12 

23  X 106 

10.39  x 106 

23.93  X 106 

Cl3 

= C12 

= Cl2 

14.94  X 106 

C23 

= C12 

= C12 

= Ci3 

C44 

6.5  X 106  = 

>4  (C11  -C12) 

6.92  x 106 

10.92  X 106 

C55 

= C44 

= C44 

= C44 

^66 

= C44 

= C44 

10.30  X 106  = 

’/>  (C11  C12) 

4 - Transversely  Isotropic 
ZincB 

5 - Transversely 
Isotropic 

6 Orthotropic 

C11 

23.35  X 106 

49.4  X 106 

34.05  X 106 

C22 

= Cli 

= Cii 

22.24  X 106 

C33 

8.85  X 106 

38.1  X 106 

21.75X  106 

C12 

4.96  X 106 

34.6  X 106 

7.05  X 106 

C13 

7.27  X 106 

9.7  X 106 

5.76  X 106 

C23 

= Ci3 

= Ci3 

5.21  X 106 

C44 

5.55  X 106 

14.2  X 106 

1.00  X 106 

C55 

= C44 

= C44 

5.00  X 106 

C66 

9.20  X 106  = 

7.4  X 106  = 

8.40  X 106 

'/i  (C11  -C12) 

'/i  (Cn  • Cl 2) 

Note:  A - Tj 

= Cjjej,  Cjj  given  in  psi 

B — Reference  1 2,  p.  278. 
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Timing  studies  of  the  interpolation  show  that  a complete  kernel  function  evaluation  requires 
.0015  cpu  sec.  using  cubic  interpolation.  Of  this  time,  .0009  sec.  is  used  for  the  interpolation 
of  the  modulation  function  and  .0006  sec.  to  calculate  U,)  and  T.-j  from  the  interpolated  data. 
The  computer  time  required  for  an  arbitrary  anisotropic  analysis  should  be  about  twice  that 
required  for  an  isotropic  analysis  of  the  same  geometry. 

VERIFICATION  OF  THE  ANISOTROPIC  BOUNDARY-INTEGRAL  EQUATION  ANALYSIS 

In  order  to  provide  final  verification  of  the  formulation  and  numerical  treatment  of  the  aniso 
tropic  three-dimensional  boundary-integral  equation  stress  analysis,  the  higher  order  program'* 
was  modified  to  use  the  interpolation  technique  for  evaluating  the  point  load  solutions.  Two 
sets  of  problems  were  run;  the  first  for  uniaxial  stress  states  and  the  second  for  stress  con 
centration  problems. 

Unixial  States  of  Stress 

Exact  solutions  exist  for  an  anisotropic  rectangular  parallelepiped  subjected  to  simple  tension 
or  pure  shear. This  allows  the  comparison  of  anisotropic  boundary-integral  equation  re- 
sults with  exact  solutions.  A single  boundary-integral  equation  map,  shown  in  Figure  1,  was 
used  for  all  of  these  test  cases.  The  map  was  rotated  to  allow  loading  to  be  applied  in  dif 
ferent  directions.  The  results  of  these  tests  are  shown  in  Tables  5 through  8. 

InTable  5 it  can  be  seen  that,  for  an  isotropic  material,  the  boundary-integral  equation  method 
gives  completely  equivalent  results  whether  closed  form  or  interpolated  kernels  are  used. 
Further,  both  sets  of  results  essentially  reproduce  the  exact  solution  to  the  problem  of  siinple 
tension.  Tables  6 and  7 show  excellent  agreement  between  exact  and  boundary-integral 
equation  solutions  for  a transversely  isotropic  and  an  orthotropic  material.  Finally,  Table  8 
gives  results  for  the  case  of  pure  shear  (t^  0,  = 0,  i = 1,  5).  The  maximum  error  for  the 

pure  shear  cases  is  less  than  .05%. 

For  the  problems  involving  uniaxial  stress  states  solutions  using  interpolated  point  load  solu- 
tions required  1.7  times  as  long  as  those  using  exact  point  load  solutions.  All  the  cases  dis- 
cussed used  linear  variation  of  boundary  data  and  the  same  precision  of  numerical  integration 
in  the  generation  of  the  equation  system. 
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TABLE  5 


Cube  in  Uniform  Tension  - (Isotropic-Material  #2) 


Exact 

BIE 

BIE 

Load 

Solution 

Numerical  Kernels 

Exact  Kernels 

Direction 

U2/Uze* 

1.00000 

1.00017 

1.00020 

Z 

Ux/U,e 

-.38983 

-.38983 

-.38983 

Uv/U,e 

-.38983 

-.38983 

-.38983 

Ux/U,e 

1.00000 

1.00018 

1.00020 

X 

Uy/Uxe 

-.38983 

-.38981 

-.38983 

^z^'^xe 

-.38983 

-.38981 

-.38983 

^y'^ye 

1.00000 

1.00018 

1.00020 

Y 

^x^^ye 

-.38983 

-.38981 

-.38983 

Uz/Uye 

-.38983 

-.38983 

-.38983 

All  displacements  are  normalized  by  the  exact  extension  for  the  applied  loading  in  Tables  5,  6 and  7. 


TABLE  6 


Cube  in  Uniform  Tension  - (Transversely  Isotropic-Material  #5) 


Exact 

Solution 

BIE 

Numerical  Kernels 

Load  Direction 

Uz/Uze 

1.00000 

1.00018 

Z 

Ux/Uze 

-.11547 

-.11554 

Uy/Uze 

-.11547 

-.11553 

1.00000 

1.00016 

X 

Uy/Uxe 

-.68464 

-.68460 

-.08028 

-.08016 

Uy/Uye 

1.00000 

1.00015 

Y 

^x^*^ye 

-.68464 

-.68462 

-.08028 

-.08016 
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Cube  in  Uniform  Tension  - (Orthotropic-Material  #6) 


Exact 

Solution 

BIE 

Numerical  Kernels 

Load  Direction 

Uz/Uze 

1.00000 

.99998 

Z 

-.12903 

-.12921 

Uy/^ze 

-.19355 

-.19365 

Ux/Uxe 

1.00000 

.99924 

X 

Uy/Uxe 

-.27000 

-.26968 

*-*z'^*^xe 

-.20000 

-.19996 

Uy/Uye 

1.00000 

1.00014 

Y 

^x^^ye 

-.17419 

-.17428 

Uz/Uye 

-.19355 

-.19357 

' TABLE  8 

Cube  in  Pure  Shear  (X-Y) 

Material  Uy/U„ 

Y 'exact 

Isotropic-numerical  .99984 

(Material  #2) 

r Isotropic-  -exact  kernels  .99993 

i (Material  #2) 

i 

r 

Transversely- - Isotropic  .99989 

(Material  #5) 

Orthotropic  1.00(X)4 

(Material  #6) 
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Stress  concentration  problems 


A second  set  of  cases  was  chosen  to  test  the  applicability  of  the  anisotropic  analysis  to  stress 
concentration  problems,  since  it  is  in  such  problems  that  the  boundary-integral  equation 
method  finds  its  engineering  application.  In  addition,  these  problems  exercise  the  quadratic 
variation  in  geometry  and  boundary  data  allowed  by  the  computer  program  used. 

The  first  problem  is  that  of  spherical  cavity  in  a cylindrical  rod.  The  geometry  and  boundary- 
integral  equation  map  used  are  shown  in  Figure  2.  Two  different  comparisons  were  made 
for  this  problem.  First  the  boundary-integral  equation  code  was  used  to  evaluate  the  stress 
concentration  for  an  isotropic  material  over  a range  of  a/w,  using  both  closed  form  and 
interpolated  point  load  solutions.  The  stress  concentration  (in  terms  of  mean  net  section 
stress)  is  compared  to  the  values  of  Reference  14  in  Figure  3.  Both  integral  equation  analyses 
show  excellent  agreement  with  Reference  14.  Of  particular  interest  is  the  fact  that  the  in- 
tegral equation  solution  using  interpolated  point  load  solutions  shows  no  more  sensitivity 

to  changing  a/w  than  that  using  closed  form  solutions.  This  indicates  that  mapping  require- 

4 

ments  for  anisotropic  analysis  should  be  no  more  severe  than  for  isotropic  analysis.  It 
should  be  noted  that,  for  this  problem,  more  precise  integration  in  the  equation  generation 

was  required  when  using  interpolated  point  load  solutions  than  when  using  closed  form  solu-  , 

tions.  This  was  also  true  in  the  notch  problem  discussed  below.  The  computer  time  ratio 
of  1.7  for  the  uniaxial  problems  increased  to  the  range  of  2 to  2.5  for  these  problems.  Some 
further  comments  on  this  question  will  be  made  in  the  last  section  of  this  paper. 

The  problem  of  the  spherical  cavity  was  also  solved  (for  a/w  = .05  ^ 0)  for  two  transversely  ■ 

isotropic  materials,  cobalt  and  zinc.  The  results  are  compared,  in  Table  9,  to  closed  form 
results  derived  for  the  limiting  case  a/w  = 0. 
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TABLE  9 


i 


Stress  Concentration  for  a Spherical  Cavity 
in  a Transversely  Isotropic  Cylindrical  Rod 


•^TN 

■^TN 

Material 

Ref.  2 

BIE-  - Numerical  Kernels 

Zinc 

1.63 

1.47 

Isotropic  {u  = .3) 

2.05 

2.03 

Cobalt 

2.31 

2.29 

Note-  - a/w  = .05  (=  0)  for  all  cases. 

As  a last  test  problem  the  stress  concentration  was  calculated  for  a deep  hyperboloidal  notch 
under  tension.  A closed  form  solution  to  this  problem  exists  for  a transversely  isotropic 
material.  The  geometry  and  boundary-integral  equation  map  for  the  problem  are  shown  in 
Figure  4.  The  results  of  the  analysis,  for  one  isotropic  and  three  transversely  isotropic  ma- 
terials are  shown  in  Table  10.  As  in  previous  cases,  the  boundary-integral  equation  solutions 
using  closed  form  and  interpolated  point  load  solutions  have  equivalent  accuracy.  Further, 
the  level  of  accuracy  is  the  same  for  both  isotropic  and  anisotropic  materials.  Finally,  it 
should  be  pointed  out  that  the  modeling  problem  in  the  notch  is  extremely  complex.  With 
the  relatively  coarse  model  used  the  results  are  sensitive  to  the  precise  location  of  the  mid 
point  nodes  indicated  in  Figure  4,  a question  which  is  independent  of  material  anisotropy. 
More  precisely,  correct  location  of  these  nodes  is  required  to  maintain  nearly  circular  cross 
sections  in  constant  z planes  and  to  ensure  that  the  BIE  model  of  the  notch  possesses  a ver 
tical  tangent  at  the  z = 0 section. 


TABLE  10 
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Stress  Concentration  in  a Hyperboloidal  Notch 


Material 

•^TN 

Ref.  3 

•^TN 

BIE  - Numerical  Kernels 

•^TN 

BIE-  - Exact  Kernels 

Zinc 

2.30 

2.42 

— 

Nickle  Base  Alloy 

3.32 

3.46 

— 

Isotropic  iv  = .3) 

3.34 

3.45 

3.22 

Cobalt 

4.02 

4.16 

Note-  - Kj|i^  = maximum  axial  stress  divided  by  average  stress  at  minimum  cross  section. 
Diameter  at  minimum  section  = 1.0,  notch  radius  of  curvature  = .1. 

It  should  be  noted  that  the  results  for  zinc  are  less  accurate  than  those  for  the  other  materials 
in  both  stress  concentration  problems.  It  is  known^®  that  certain  mathematical  properties  of 
the  point  load  solution  for  zinc  differ  from  those  of  the  other  materials  used  in  this  study.  It 
is  believed  that  the  effect  of  this  difference  in  the  present  formulation  is  to  require  more 
accurate  numerical  integration  of  (4)  for  zinc  to  achieve  the  same  accuracy  in  the  stress  con 
centration  factors. 
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CONCLUSIONS 

It  has  been  shown  that  anisotropic  point  load  solutions  can  be  numerically  evaluated  with  es- 
sentially arbitrary  accuracy.  Further,  accurate  boundary-integral  equation  stress  analysis  can 
be  carried  out  using  these  solutions.  Use  of  an  appropriate  interpolation  technique  makes  the 
anisotropic  analysis  feasible  for  use  in  engineering  applications. 

Future  work  could  profitably  be  directed  at  further  increases  in  computing  efficiency. 

In  particular: 

1.  It  was  noted  in  solving  the  stress  concentration  problems  that  more  accurate  integration 
was  required  when  using  interpolated  point  load  solutions.  It  is  likely  that  this  is  re- 
lated to  slight  oscillations  in  the  modulation  function  induced  by  the  interpolation 
scheme.  Further  investigation  should  allow  at  least  partial  resolution  of  the  problem, 

with  the  possibility  of  a significant  saving  in  solution  time. 

2.  There  may  well  exist  a better  means  of  reading  the  modulation  function  data  base. 

The  presently  used  Lagrange  interpolation  is  a standard  technique  which  incorporates 
no  knowledge  of  the  physics  of  the  problem. 

3.  A closed  form  point  load  solution  exists  for  transversely  isotropic  materials^®.  For 
this  material  class  use  of  the  closed  form  solution  may  lead  to  significant  time  savings. 

Finally,  more  work  remains  to  be  done  in  extending  the  anisotropic  boundary-integral  equa- 
tion analysis  to  include  thermal  and  rotational  body  forces,  both  of  substantial  importance 
' in  engineering  applications. 
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FIGURE  1 


BOUNDARY-INTEGRAL  EQUATION  MAP  FOR 
UNIAXIAL  STRESS  STATE  PROBLEMS 
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FIGURE  2 BOUNDARY-INTEGRAL  EQUATION  MAP  FOR  SPHERICAL  CAVITY 


FIGURE  3 STRESS  CONCENTRATION  FOR  A SPHERICAL  CAVITY 


